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We perform an analytical sensitivity analysis for a model of a continuous-time branching process evolving 
on a fixed network. This allows us to determine the relative importance of the model parameters to the 
growth of the population on the network. We then apply our results to the early stages of an influenza¬ 
like epidemic spreading among a set of cities connected by air routes in the United States. We also 
consider vaccination and analyze the sensitivity of the total size of the epidemic with respect to the 
fraction of vaccinated people. Our analysis shows that the epidemic growth is more sensitive with respect 
to transmission rates within cities than travel rates between cities. More generally, we highlight the fact 
that branching processes offer a powerful stochastic modeling tool with analytical formulas for sensitivity 
which are easy to use in practice. 
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1. Introduction 

Branching processes are powerful stochastic models that describe the evolution of populations of individ¬ 
uals which reproduce and die independently of each other according to specific probability laws. They 
are playing an increasingly important role in models of population biology including molecular biology, 
ecology, epidemiology, and evolutionary theory [6], as well as in other scientific areas such as particle 
physics, chemistry, and computer science [21]. Typical performance measures of these models include the 
distribution of the instantaneous and cumulative population sizes at a given time, the extinction probability 
of a population, and its asymptotic growth rate and composition. 

In order to apply branching processes to real-world problems, the parameters of the model must be 
estimated from available data sets. Small errors or changes in the parameters may lead to notably different 
model outputs. A sensitivity analysis may quantify the impact of each parameter on the performance 
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measures of the model. This may be useful if we want to know which parameters influence the growth of the 
branching process the most. To the best of our knowledge, sensitivity analysis has received scant attention 
in the branching processes literature. Apart from a few papers dealing with very specific questions such as 
[14] and [16], there is no complete study of the topic for more general classes of branching processes. 

Here we consider a continuous-time Markovian branching process evolving on a fixed network. Such 
a process belongs to the class of Markovian trees, which can be seen as particular multitype branching 
processes, and offer considerable modelling versatility and appealing computational properties, see [12], 
[23] and [27]. The concept of a branching process evolving on a network is certainly not new, but our aim 
in this paper is to present branching processes and their properties in an accessible way, and to show how to 
derive practical, analytic sensitivity formulas for the main quantities of interest. For that purpose, we first 
discuss the typical performance measures of the model (as listed earlier), and we then perform an analytical 
sensitivity analysis of each performance measure with respect to the model parameters. As opposed to the 
performance measures themselves, most of the sensitivity formulas have an explicit expression that can be 
easily used in practical situations. 

We then illustrate our theoretical sensitivity results on a topical real-world problem: the spread of an 
influenza-like epidemic on a network of cities in the United States through air traffic. It is well known 
that the early stages of an epidemic can be approximated by a branching process [10], [35]. The averaged 
branching process is essentially the linearisation of the nonlinear Susceptible-Infected-Removed (SIR) com- 
partmental model around the disease-free equilibrium. We refer to Hethcote [28] and the book of Keeling 
and Rohani [30] for a good introduction to the modeling of infectious diseases, to Arino [5] for an overview 
on diseases in metapopulations, and to Balcan et al. [9] for a recent computational model for the spatial 
spread of infectious diseases. 

The sources of errors in the outputs of a model of epidemic evolving on a network are manifold. They 
may arise either from the simplifying assumptions of the model, which necessarily neglect many features 
of the real world, or from imprecise epidemiological or mobility data. The former is addressed in the most 
sophisticated models by taking into account a large number of explanatory variables or compartments, such 
as age, gender, location, etc., and modeling the dynamics in a complex nonlinear way. The mobility data, 
at a global level or at the level of a large country such as the United States, is often estimated from the 
number of passengers flying from one airport to the other. These data are easily available and are assumed 
to account for a large part of the mobility. A discussion on the relevance of these data is included in [8]. The 
epidemiological data are however much harder to estimate accurately. For instance, the average number of 
infectious contacts that an individual infected by seasonal influenza makes in a day have been estimated to 
range from 0.55 to 1.44, see [18]. 

Sensitivity of epidemic models to parameters has already been studied for different types of diseases: 
Hyman and LaForce [29] studied the sensitivity of a multi-city deterministic epidemic model for influenza 
with respect to the epidemiological data. They assumed that the transmission and recovery rates are the 
same in each city and show that the recovery rate is the most important single parameter. Chitnis et al. [15] 
perform a sensitivity analysis on a non-linear compartmental model of malaria transmission, and obtain an 
analytical formula for the basic reproduction number Rq and for its sensitivity. 

In this paper, we compute the sensitivity of the branching process approximation with respect to both 
epidemiological and domestic mobility data. We find that the results are significantly more sensitive with 
respect to the epidemiological data than with respect to the domestic mobility data. This provides confi¬ 
dence in the use of the general approximation of domestic mobility by airport traffic, but suggests that the 
most stringent limitation of current models is in the relatively imprecise epidemiological data, rather than 
in the sophistication of the model or the estimation of domestic traffic. This is in contrast with the role of 
international air travel, which can have a more important impact on the development of some diseases such 
asHlNl [31]. 

Finally, we refine our sensitivity analysis with respect to the epidemiological parameters by studying 
the effect of a vaccination campaign. The sensitivity analysis of the epidemic size with respect to the 
vaccination rate allows us to compute a practically relevant quantity: the impact of each supplementary 
vaccination shot on the total number of infected people, when the vaccination has succeeded in stopping 
the exponential spreading of the disease. In this regime the stochastic branching process approximation is 
valid at all times [10]. 

In an appendix, we also introduce the possibility of on-board transmission. Infection during transporta- 
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tion has been analyzed from a theoretical point of view for the SIR model in [32], but the authors have 
not applied the model to real situations. Studies on the transmission rates on board airplanes exist (see for 
instance the report of the European Centre for Disease Prevention and Control [1]), however, to the best of 
our knowledge, data are not sufficient to provide accurate on-board transmission functions for influenza. 
Therefore we assume that the transmission probabilities are given by a specific function of the flight time, 
and we show how this assumption affects the sensitivity results. We show in particular that considering 
on-board transmission can noticeably modify the shape of the sensitivity curves. 

The paper is organised as follows. In Section 2, we describe the model of branching process evolving 
on a network and we discuss several performance measures which can be computed from the model. In 
Section 3, we define the sensitivity and the elasticity of a model output with respect to a parameter, and 
we derive analytical expressions for the sensitivity of each performance measure discussed in Section 2. 
We end the section by an illustration of the sensitivity results on a simple artificial example. In Section 4, 
we apply the sensitivity results on a model for the early spread of an influenza-like epidemic on a network 
of cities in the United States. Finally, we introduce vaccination and perform a sensitivity analysis with 
respect to the proportion of vaccinated people. We conclude our paper in Section 5. The technical details 
and supplementary information are provided in six appendices. 

Throughout the paper, column vectors are denoted by x and row vectors are denoted by The 
column vector e, denotes the unit vector with a 1 at the /th entry and zeros elsewhere, and 1 and 0 are 
column vectors of which all elements are respectively equal to one and zero; the size of these vectors is 
generally clear from the context. 


2. A model of branching process evolving on a network 

We consider a fixed network represented by an undirected graph where = {1,2,..., n} is the set 

of nodes and S is the set of vertices. Two nodes i and j are adjacent when either i = j or the edge [i,]) 
belongs to 

We define a continuous-time stochastic model of a population evolving on the network as follows. The 
process starts at time 0 with a single individual located at node i € One of the following three types of 
event may then happen to this individual while at node i: 

• The individual moves from i to an adjacent node J ^ i; this happens at rate (that is, Tij transitions 
from i to J occur on average per time unit, per individual at node i). 

• The individual gives birth to k children (k ^ 1) and simultaneously moves to node j adjacent to i, 
the k children respectively starting their life at nodes ji, j 2 , ■ ■ ■ ,jk adjacent to /; this happens at rate 

{Bk) 

i'Jl jl-JkJ' 

• The individual dies; this happens at rate di. 

All individuals living on the network behave independently of each other with the same rules as the initial 
individual. The transition rates, birth rates and death rates are gathered respectively in a n x n matrix 
T — {Tij), a sequence of n x matrices - Aj)’ ^ an n x 1 vector d ~ {di). The 

diagonal elements of T are strictly negative and |7],j is the parameter of the exponential distribution of 
the sojourn time of an individual at node i before one of the three abovementioned events occurs. These 
elements are computed such that the matrices and vector satisfy T\ +<7 = 0. The resulting 

Markovian population process is a branching process characterized by the set of matrices {T 
belonging to the class of Markovian trees, see [12] and [26]. 

Example 2.1 We illustrate the structure of the matrices T, d on the simplest network with 

n = 2 (adjacent) nodes. The 2x2 transition rate matrix T is then given by 


Til Ti2 
721 722 


where the diagonal entries 7], will be described further. Assume that the individuals of the branching 
process can give birth to at most two children at each birth event. This means that the birth rate matrices 
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Bk are nonzero only for A: = 1 and k = 2. The 2x4 matrix Bi has the following structure 


Bi = 


(fil)l,22 

(^1)2,11 (^1)2,12 (51)2.21 (51)2.22 


and the 2x8 matrix B 2 has the following structure: 


52 = 


(52)1.111 (52)1,112 (52)1,121 (52)1,122 
(52)2,111 (52)2,112 (52)2,121 (52)2,122 


(52)1.211 (52)1,212 (52)1.221 (52)1,222 
(52)2.211 (52)2,212 (52)2.221 (52)2,222 


For instance, the entry ( 52)1 122 is the rate at which a parent at node 1 gives birth to two children and 
instantaneously moves to node 2, while one of his children stays at node 1 and the second one moves to 
node 2. Note that the indices in the entries of Bi and B 2 are ordered lexicographically by convention. 
Finally the 2x1 death rate vector is af = [d\^d 2 \^. In order for T\ +B 1 I +521 +d = 0 to hold, the 
diagonal entries of the matrix T are then given by 


511 


T 22 


—d\ — Ti2 


—c /2 — 521 


2 


2 


H (^i)i,;ii ^ (52)i,;ij2;) 


2 


2 


52 ( 51 ) 2 ,;i2 52 (^ 2 ) 2 ,;, jj;- 


In the next subsections, we describe several performance measures that can be computed for a branching 
process evolving on a network. For the clarity of the presentation, the performance measures will be 
discussed first in the scalar case (n = 1) which corresponds to the standard Markov branching process [22, 
Chapter V], and next in the matrix case for an arbitrary number of nodes n > 1. Note that in the scalar 
case, the matrices T, {5i^}i.^i and d become all scalar: Bj^ is the rate at which an individual gives birth to k 
children, d is the death rate, and T = —d — Y^k^k- 

The matrix case will require the use of the Kronecker product between matrices, which is defined as 
follows: for an n x m matrix A and a px q matrix B, the Kronecker product of A and B, denoted by A 0 B, 
is an np X mq matrix defined as 


Aii5 

Ai25 . 


A2i5 

^225 . 


A„i5 

An2B . 



It will also require the knowledge of the Perron-Frobenius Theorem for nonnegative matrices, see [37]. 

Most of the results in the scalar case can be found in [22], and most of the results in the matrix case can 
be found in [23] and [26] for the case where Bk = Q for k^2 (the Markovian binary tree case). 


2.1 Instantaneous population size 


Scalar case n = 1. Let Z(f) denote the population size in the branching process at time t. The process 
{Z(f),f ^ 0} is a continuous-time Markov chain on the nonnegative integers, where state 0 is absorbing 
and all other states are transient (that is, they will be visited a finite number of times with probability one). 
The so-called generating function of Z(f) is defined as Q{t,s) = Lyt>o5[Z(f) = k]s^, where |s| ^ 1, and is 
well known to satisfy the backward Kolmogorov differential equation 


BQ{t,s) 

dt 


d + TQ{t,s)+Y.BkQ{t,sf+\ 

k^\ 


( 2 . 1 ) 


with 2(0, s) = s if the population starts with a single individual at time 0 (which we assume here). 

All the moments of Z(f) can be obtained by successive derivatives of Q{t,s) at s = 1. In particular, the 
mean population size at time t, M{t) = E[Z(t)], is given by M{t) = ((92(f,s)/(9s)|.v=i. By differentiating 
(2.1) with respect to s, at s = 1, we obtain the linear differential equation for M{t) 


dM{t) 

dt 


Q.M{t), 


( 2 . 2 ) 
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withM(O) = 1, where 

= r + ^ Bj^ (fe + 1) = —d + ^ kBh^. 
k^\ k^l 

By solving (2.2) we obtain that the mean population size at time t is given by 


M{t) = exp(f2f). 


(2.3) 


(2.4) 


Matrix case n > 1 . We now extend the above definitions and results to the n-dimensional case. The 
vector Z(f ) = [Zi (f ),..., Z„ (f)] records the population size at time t at each of the n nodes. The evolution of 
the branching process now depends on the node occupied by the initial individual at time 0. We therefore 
define the conditional probability generating function of Z(f), given that the process starts with a first 
individual at node i, as 

k^O 

where k = {ki,... ,k„) € N", s = (si,... ,s«), |s,| ^ 1, and s* := s** • • As shown in [26], the vector 
function Q(t,s) = (Qi(t,s)) satisfies the matrix analogue of (2.1), 


dQ{t,s) 

dt 


d + TQ{t,s)+Y,BkQ{t,sf+^'>, 

k 


(2.5) 


with 0(0,s) = s, where stands for the {k-\- l)st-fold Kronecker product of the vector Q{t,s) 

with itself; = 1, and 0 Q{t,s), for k^\. 

The mean population size at time t is now given by a matrixM(f) = {Mij{t)), where My(f) =£[Z,(f)|Z(0) 
e,] is the (conditional) mean number of individuals at node j in the population at time f, given that the 
process started at time 0 with one individual at node i. The entries of M{t) being obtained as = 

{dQi{t^s)/dsj)\s=\, the same differential equation (2.2) and solution (2.4) hold for the matrix M{t), the 
only difference being that M(0) = /, and Q is now a matrix given by 



Note that the matrix exponential is defined as 


e 


nt 


E 


n\ 


( 2 . 6 ) 


2.2 Cumulative population size 

Scalar case n = 1. Let N{t) be the cumulative number of individuals born in the branching process 
until time t, and let G{t,s) be the probability generating function of N{t). Similar to Q{t,s), G{t,s) satisfies 
the differential equation 


dG{t,s) 

dt 


ds + TG{t,s)+Y^BkG{t,sf+\ 
k 


(2.7) 


with G(0,s) = s. Note that the only difference with Eq. (2.1) is that here d is multiplied by s. 

Again, all the moments of N{t) can be obtained by successive derivatives of G{t,s) at s = 1, and 
we focus here on the mean cumulative number of individuals born in the branching process until time f, 
D{t) = £[A(f)] = ((9G(f,s)/(9s)|s=i. By differentiating Eq. (2.7) with respect to s, at s = 1, we obtain the 
following differential equation for D(t): 


dD{t) 

dt 


d + Q.D{t), 


with Z)(0) = 1 and where the scalar Q is defined in (2.3). The mean cumulative population size until time 
t is thus given by 


D(f) = [/ —exp(f2f)] (—f2) * c/ + exp(f2f). 


(2.8) 
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Matrix case n > 1. In the n-dimensional case, we extend the definition of G{t,s) to its vector analogue; 
we define the conditional probability generating function of N{t), given that the process started with a first 
individual at node i, by 

Giit,s) = Y,P[N{t)=k\Z{0)=ei]s‘‘. 

Note that here, s is still scalar. The matrix analogue of Eq. (2.7) for the vector G{t,s) = {Gi{t,s)) is 

= ds + TGit,s)+Y,BkGit,sf+^\ (2.9) 

k 

with G{Q,s) = s 1, see [26]. Since we now condition on the node of the initial individual at time 0, the mean 
cumulative number of individuals born in the population until time t becomes a vector D{t) = {Di{t)) where 
Di{t) = £'[N(f)|Z(0) = e,] = ((9G,(f,s)/(9s)|.v=i. Using the same technique as in the scalar case, we find 
that the vector D{t) is given by the matrix analogue of (2.8), 

D{t) = [I-exp{Qt)]{-Q)-U + exp{Qt)l, (2.10) 

where the matrix Q is defined in (2.6). Note that this expression requires Q to be nonsingular, which is 
generally the case in most practical situations. 


2.3 Extinction probability 

Scalar case n = 1. Due to the transient nature of the strictly positive states of {Z(f)}, the population 
in the branching process either eventually grows without bound or becomes extinct, there is no stationary 
behavior (but for trivial cases). We denote by q the probability of eventual extinction of the process, that 
is, q = lim;^ooP[Z(f) = 0]. Since P[Z{t) = 0] = Q{t,Q), an equation for q can be obtained by taking s = 0 
and f —> oo in Eq. (2.1): 


0 = d + Tq + Y.Bkq^""'^- ( 2 . 11 ) 

k 

This non-linear equation has potentially more than one solution, but it can be shown that q is its minimal 
nonnegative solution. 

Matrix case n > 1. The n-dimensional version of q is the vector q = {qi) where qi = limr^„o7’[Z(f) = 
0|Z(0) = e,] is the conditional extinction probability of the branching process, given that it starts with one 
individual at node i. Similar to the scalar case, we can show that the vector q is the (componentwise) 
minimal nonnegative solution of the matrix analogue of (2.11), 

0 = d + Tq + Y.Bkq^''^'^. ( 2 . 12 ) 

k 

Several algorithms with a probabilistic interpretation have been developed to solve for q, see [12], [25], 
and [27]. 

2.4 Extinction criteria 

Scalar case n = 1. Various criteria can be used to determine whether the population eventually 
becomes extinct with probability one or has a positive probability of growing without bound, that is, 
whether ^ = 1 or ^ < 1. We shall consider two equivalent extinction criteria here. 

Eirst observe that, by (2.4), the mean population size M{t) has different asymptotic behavior depending 
on the sign of Q defined in (2.3): as f —> oo, M(t) —> 0 if <0 (subcritical case), M(t) = 1 if = 0 
(critical case), orM(f) —)• oo if >0 (supercritical case). It follows that the population eventually becomes 
extinct with probability one {q = 1) if and only if ^ 0. 

Another threshold quantity is given by the mean number of children generated by an individual during 
its entire lifetime, that we denote by R. The population has a positive probability of surviving {q < 1), if 
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and only if R> 1, that is, if on average, every individual is replaced by more than one individual. In the 
scalar case, R is the ratio of the total rate at which an individual generates a child to the death rate, that is, 

R ^ LkkBk 
d 

Note that the two criteria are indeed equivalent since Q. ^ 0 ^ 1. 

Matrix case n > 1. The two extinction criteria described in the scalar case have a counterpart in the 
n-dimensional case. Since Q is now a matrix (as defined in (2.6)), the condition ^ 0 for almost sure 
extinction needs to be adapted; the eigenvalue Qq of maximal real part of Q (also called the Perron- 
Frobenius eigenvalue) now plays the role of the threshold quantity: 

Q = f ^ ^ 0 . 

We define the matrix R = (Rij), where the entry Rij is the expected total number of children born at 
node j from a parent born at node i, during the entire lifetime of the parent. The second extinction criterion 
then relies on the eigenvalue Rq of maximal real part of the matrix R: 

q=l Rq^I. 

An explicit expression for R, which generalizes (2.13) to the matrix case, is given in Proposition A.l in 
Appendix A. 

The matrix R and its dominant eigenvalue Rq will play a fundamental role in the epidemic application, 
as shown in Section 4. 


2.5 Asymptotic node frequency 

It may be interesting to know what proportion of the inviduals living in the network at time t are located 
at node i. We also call this proportion the frequency of node i at time t, and when t —)■ oo we talk about 
the asymptotic node frequency. Note that these notions make sense only when the network has at least two 
nodes; therefore, we only consider the matrix case in this section. 

The matrix exponential exp(f2) has the Perron-Frobenius eigenvalue exp(f2o). Hence, as a conse¬ 
quence of (2.4) and Perron-Frobenius theory, the expected population size at time t is asymptotically given 
by M{t) ~ exp(f2of) as f —> oo, where and v are respectively the left and right Perron-Frobenius 
eigenvectors of exp(f2); this holds provided that the dominant eigenvalue has multiplicity one, which is 
generally the case in most practical situations. This implies that for all i,j, 


Mij{t) uj 
{M{t)\)i ^ 


as f —> oo. 


(2.14) 


and we may use (2.14) to obtain the proportion of the population living at each node of the network as time 
goes to infinity. 


3. Sensitivity analysis 

In practical situations, precise values for the parameters of the model can be difficult to 

obtain. It is therefore important to perform a sensitivity analysis of the model with respect to perturbations 
or errors in the parameters. This section addresses the sensitivity analysis of the measures of interest 
defined in the previous section with respect to the parameters of the model. 

Let p be a parameter of the model (for instance p = 7)/, p = or p = di), and let X be a measure 

of the model (for instance X = Mij{t), X = Dfr), or X = Rq). The sensitivity of X with respect to p is 
defined by the local slope of X, considered as a function of p: 
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The scale of X and p may be different; it is therefore convenient to consider proportional perturbations 
instead of absolute ones. The proportional response to a proportional perturbation is the elasticity (also 
called sensitivity index in the context of mathematical epidemiology [15]). The elasticity of X with respect 
to p, denoted by , is defined by the ratio of the relative increase of X to the relative increase of p: 


d log p ^ X 


(3.2) 


The interpretation of the elasticity is as follows: if 1^ = a, it follows that if p increases by r%, then X 
increases (or decreases, depending on the sign of a) by approximately 


100[exp(fllog(l+r*0.01))-l]% (3.3) 

(see Appendix B for more details). Note that for r = 1, log(l.Ol) sa 0.01 so that 100[exp(fl/100) — 1] « a 
when \a\ is small. 

We now derive analytical formulas for the sensitivity of the performance measures defined in the pre¬ 
vious section. Here we focus on the matrix case only, as this is the case leading to relevant discussions. 


3.1 Sensitivity of the instantaneous population size 

In order to characterize the sensitivity of the population size generating function Q{t,s), we differentiate 
(2.5) with respect to p. This leads to a matrix linear differential equation for dpQ{t,s) which does not have 
any closed-form solution, see Eq. (A.3) in Appendix A. 

Recall from (2.4) that the mean population size matrix is M{t) = exp(f2f). Derivatives of the matrix 
exponential have been investigated by several authors, see for instance [4], [33] and [34]. We present here 
a simple method which provides an exact expression for dpM{t) and involves the computation of a matrix 
exponential of size 2n x 2n. Let us consider the system of differential equations 

r dM{t) = 

\ dtdpMf) = QdpM{t)+dpQM{t), ^ ’ 

where the first equation is Eq. (2.2) satisfied by M(f), and the second equation is obtained by differentiating 
the first equation with respect to p. The system (3.4) may be equivalently rewritten as 


■ dpM{t) ' 


■ Q. 

dpQ, 


■ dpM{t) ' 

M{t) 


0 

Q. 


M{t) 


with initial condition [dpM{0)^M{0)Y = [0,/]^, where 0 and I denote the n x n zero matrix and identity 
matrix respectively. This is a new differential equation of the form dtY (f) = AY{t), of which the solution 
is given by Y (f) = exp(Af)y(0). Therefore, 


where 


dpM{t) = [7,0] - exp 


Q. 

0 


dpQ, 

Q. 


0 

I ’ 


dpQ=dpT+Y. dpBk 0/0 1 

jt^l i=0 


(3.5) 


(3.6) 


A second method to compute dpM{t), which has the advantage of not requiring any matrix exponential 
computation, is provided in Appendix A. 


3.2 Sensitivity of the cumulative mean population size 

Eocusing now on the sensitivity analysis of the cumulative population size, we observe that, similar to 
dpQ{t,s), dpG{t,s) satisfies a linear matrix differential equation with non-constant coefficients which does 
not have any closed-form solution. 
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The explicit formula for dpD{t) follows easily from the one for dpM{t). Recall the expression for D{t) 
given in (2.10) and note that 

dp{-a)-^ = {-a)-^dpa{-a)-\ 


We therefore obtain 

dpD{t) = -dpM{t){-Q.)-'^d+[I-&x^{Q.t)]{-Q.)-'^dpQ.{-Q.)-'^d + dpM{t)\, (3.7) 
where dpM{t) is computed in Section 3.1. 


3.3 Sensitivity of the extinction probability 

Here we derive an explicit expression for dpq in terms of the extinction probability vector q. By differen¬ 
tiating (2.12) with respect to p, we obtain 

0 = dpd + dpTq + Tdpq + Y^dpBtq^’^^''^ 

k k i=0 


SO that, in the non-critical case, 


dpq = -d> ' 



(3.8) 


where 

k 

0 = r + ^). 

k 1=0 

Note that when q= 1, 0 = £2. The non-singularity of the matrix 0 in the non-critical case is ensured by 
the next proposition. 

Proposition 3.1 If the Markovian tree is supercritical or subcritical, then the eigenvalues of 0 all have a 
strictly negative real part. In the critical case, 0 is singular. 

We omit the proof which follows the same lines as in [26, Theorem 6]. 


3.4 Sensitivity of the extinction criteria 


Recall that the eigenvalues of maximal real part of Q. and R, denoted by and Rq respectively, are 
key quantities to determine whether the branching process almost surely becomes extinct or not. Even 
though no explicit expression can be written for £2o and Rq, analytical expressions can be derived for their 
sensitivities dpQo and dpRo, as we show now. 

Let A be any generic nx n matrix (that is, A represents both £2 and R). Let Aq be the eigenvalue of 
maximal real part of A, and let and v be the left and right eigenvectors of A corresponding to Aq, scaled 
such that u^v = 1. Then, 


n n 


dpAo = Y^Y^ 
<= 1.7=1 


BAq dAjj 
dAij dp ’ 


where dAo/dAij = UiVj [13, Chapter 9], so that 


dpAo = dpAv. 


Linally, when A = £2, dp£2 is given in (3.6), and when A = R, we use the explicit expression (A.l) for R to 
obtain an expression for its derivative dpR, see (A.2) in Appendix A. 
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Eig. 1. Graph of a network with n = 4 nodes. 


3.5 Sensitivity of the asymptotic node frequency 

Recall from Section (2.5) that and v are the Perron-Frobenius left and right eigenvectors of exp(f2) 
associated with the dominant eigenvalue exp(f2o). Let {{Xi,uJ ,Vi)^ 1 ^ ^ n} be the full set of eigenvalues 
of exp(f 2 ), with their corresponding left and right eigenvectors, with Xi — exp(f 2 o), uj ~ and vi = v, 
and assume that the pairs of eigenvectors are normalized such that uJ v, = 1 for all i. Then, the sensitivity 
of the asymptotic node frequency is given by 


mTI 


dpU~^ (u~^ 1) - U~^dpU~^ 1 


where 


j " u'{dpexpQ)vi j 

dpU' = 1 ^ -^ U , 

hi exp(f2o)-A,- 


(3.9) 


(3.10) 


which is obtained following the same arguments as in [13, Chapter 9], and where dpSxpQ is computed in 
Section 3.1. 


3.6 Illustrative example 

In this section, we illustrate the calculation of a few sensitivity results on a branching process evolving on 
the simple network with n = 4 nodes depicted in Figure 1. Assume that every individual generates a single 
child at each birth event, that is, 8 /^ = 0 for kf^2.To simplify the notation, we drop the index 1 in the birth 
rate matrix Bi. The structure of the network induces the following structure for the matrices and vector 
characterizing the branching process: 



■ T’li 

Tn 

Tn 

7’i4 


d\ 

T = 

721 

T^ii 

T 22 

0 

0 

733 

0 

0 


d 2 

d^ 


741 

0 

0 

744 


d4 


Byn 

71i;12 

71i;13 

71i;14 

71i;21 

71i;22 

71i;23 

71i;24 

71i;31 

^l;32 

71i;33 

71i;34 

71i;41 

71i;42 

71i;43 

71i;44 

82-11 

7l2;12 

0 

0 

7l2;21 

B 2;22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Bi-yi 

0 

7l3;13 

0 

0 

0 

0 

0 

^3:31 

0 

7l3;33 

0 

0 

0 

0 

0 

_7l4;ll 

0 

0 

7l4;14 

0 

0 

0 

0 

0 

0 

0 

0 

7l4;41 

0 

0 

Z?4;44 _ 


where the diagonal of T is such that Tl + Bl + d = 0. To further simplify the model, we set di = d 
for all i, B\-i j — bi for all i,j, B 2 jj = ^2 for 67 = 1)2, Byjj = ^3 for i,j = 1,3, and B 4 ;,./ = ^4 for 
i,j = 1,4, for some parameters ( 7 , 7 >i, 7 > 2 ,^ 3 ,^ 4 . In this special case, the expression (2.6) for Q reduces to 
Q = T +B(/C) 1 + I (g)/) where 


T = 


-Ti2 — 7 ’ i 3 — Ti4 — d— I6bi 
721 
T 31 
741 


Tn 

-T 21 —d — 4 b 2 
0 
0 


Tn 

0 

—731 — d — 4^3 

0 


7’i4 

0 

0 

-Till —d — 4b4 
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Table 1. Sensitivity analysis of the toy example for a given set of parameters. 


p 

Value of p 

r"n{2) 

rOi(2) 

-yRo 


Tn 

2 

-0.4683 

-0.3764 

-0.0175 

0.0198 

Tn 

6 

-0.4679 

-0.3126 

-0.0133 

0.0151 

Tu 

8 

0.6255 

0.6723 

0.0307 

-0.0298 

Tn 

1 

0.1564 

0.1247 

0.0065 

-0.0062 

Tn 

1 

0.0733 

0.0484 

0.0023 

- 0.0022 

Tn 

1 

-0.0826 

-0.0891 

-0.0044 

0.0034 

d 

10 

- 20.0000 

-17.9012 

-0.9804 

1.0131 

b\ 

1 

7.4966 

6.9555 

0.2608 

-0.3174 

b 2 

2 

3.1233 

2.7226 

0.1227 

-0.1286 

bi 

3 

6.5588 

5.7174 

0.2448 

-0.2451 

b4 

4 

9.9944 

8.7123 

0.3478 

-0.3220 


and 


B(/(g)l + l(g)/) = 


8^1 

8bi 

8^1 

8^1 

4^2 

4b2 

0 

0 

4^73 

0 

4^3 

0 

4b4 

0 

0 

4b4 


We assume that the branching process starts with one individual at node 1 at time f = 0, and we are 
interested in the population size at time t = 2 and in the extinction probability. With the specific parameter 
values provided in the second column of Table 1, we obtain 


Mil (2) =308.7914, Di (2) = 4468.1, /?o = 1-3426, ^i =0.7390, 


so the process is supercritical. The elasticities of the above four quantities with respect to each model 
parameter are shown in Table 1. They can be interpreted using (3.3): for instance, if we increase the value 
of ri 4 by 10%, then the value of Mu (2) increases by approximately 6.14% and the value of qi decreases 
by approximately 0.28 %. We see that the parameters which influence the most the population growth are, 
in decreasing order: d, ^ 4 , bi, b^, b 2 , 7 ) 4 , 7 ’i 2 , 7 ) 3 , 72 i, T 41 and 731 . We also see that some transitions 
between the nodes positively influence the population growth (e.g. from 1 to 4), while other transitions 
have a negative influence (e.g. from 1 to 2 ). 


4. Application: Sensitivity analysis of an inlluenza-like epidemic spreading on a network of cities 

In this section we apply the sensitivity results developed in the previous section to a stochastic model for 
the initial spread of an influenza-like epidemic among a set of cities connected by air routes in the United 
States. This allows us to determine the relative importance of the model parameters to the disease spread. 
For this purpose, the early stages of the epidemic are approximated by a branching process evolving on a 
network as described in Section 2. 

Appropriate branching processes are frequently used to approximate the process of infection during the 
early stages of a general stochastic epidemic in a large closed and homogeneously mixing population. In 
theory, during the course of a major epidemic, the epidemic grows like a branching process until about \/N 
members of the population of N individuals become infected, see [11]. The accuracy of the approximation 
actually depends on our own error criteria and is discussed in Appendix C. As shown in Section 3, the 
tractability of branching processes make them perfectly suitable for a sensitivity analysis. 

4.1 Description of the model 

We assume a constant, large, and homogeneously mixing population in each city. We use a branching 
process on a network, characterized by a set of matrices { T, i,d} that we detail below, to model the 

spread of the disease among the cities during the early stages of the epidemic. As discussed in Appendix C, 
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our model of branching process is suitable as long as the number of susceptible individuals is much larger 
than the number of infected individuals in each city, which justifies why we focus on the early stages of the 
disease. In the rest of this section, the time unit is the day. 

Individuals in the branching process correspond to people infected by the disease, and each of the n 
nodes of the network corresponds to a city in the U.S. An individual makes a transition from node i to node 
j in the network when he travels by plane from city i to city j. Cities are connected through a symmetric 
n X n air travel matrix A where the entry A,y is the average number of passengers per day from city i to city 
j for i 7 ^ i (by convention we set A„ = 0). The symmetry assumption for A implies that the population of 
each city remains constant, see for instance [15] and [18]. The travel rate per individual per day from city 
i to city j is obtained by dividing the daily average number of passengers going from city i to city j by 
the population of city i, that we denote by A, ; the nondiagonal entries of the transition rate matrix T are 
therefore given by 


T,.i = i¥^j- (4.1) 

For the purpose of our application, we call the matrix T the travel rate matrix. 

In a first approach, we assume that infected people transmit the disease to new individuals within cities 
only, not during their plane travel. An individual at node i infects a new individual at rate j3,; this parameter 
thus corresponds to the average number of infectious contacts per day by an infected individual in city 
i. The rates j3, depend on the cities and are gathered in a transmission rate vector ]3. Since only one 
transmission may occur at a time, the sequence of birth rate matrices in the branching process is 

such that the only nonzero entries of the matrix B\ are 


(Bi )(■;«= A for 


and 

Bi^ = 0 fox 2 . 

Finally, infected individuals may be “removed” from the population by recovering (or dying) at rate di 
in city i, which corresponds to the inverse of the mean time of contagion of an infected individual in city 
i. These rates correspond to the death rates in the branching process and form the removal rate vector d. 
Recall from Section 2 that the diagonal elements of the matrix T are then such that 7’l+Bil+«/ = 0, that 
is, Tii = -pi - di - Y.j^i Tij. 

Due to the simple structure of the matrix Bi, we have Bi(l(g)/) = Bi(/(g)l) = diag(/3), so that the 
expressions for the matrices £2 and R, given in (2.6) and (A.l) respectively, simplify to 

f2 = r + 2diag(j3), B = -[/ + r-'diag(/3)]^‘ T-'diag(/3). 

In the epidemic context, R is the matrix of mean number of secondary infections generated in each city 
by an average infected individual in an entirely susceptible population. The dominant eigenvalue Rq of R 
is called the basic reproduction number, which is a key quantity for determining whether an infection can 
invade and persist with a positive probability: the infection can invade the population if and only if Bq > 1 
(supercritical case). When Bq < 1 (subcritical case) the disease simply dies out, and when Bq = I (critical 
case) the disease becomes endemic, meaning that the proportion of infected individuals remains constant 
over time. 

In a second approach, we take into account possible transmission on board airplanes. The description 
of the model in this case is given in Appendix D.l. 

4.2 Data 

Air travel data on the daily average number of passengers for each city-pair {i, j) among a group of n = 114 
cities are obtained from the Domestic Airline Fares Consumer Report of the US Department of Transporta¬ 
tion for the first Quarter of 2011 [2]. This report provides information on the 1,000 largest city-pair markets 
in the 48 contiguous states, among which the number of one-way passenger trips per day. These markets 
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Table 2. Travel rates from city i (which we assume here to be New York) to city j (Tij, travel frequency per individual per day from i to 
j), transmission rates in city j (fij, average number of new infection per individual per day in city j), and mean cumulative epidemic 
size after 14 days given that the disease started in city j (Dy(14)), for four selected cities j. 


City j 

T,i 


D,(14) 

New York 

- 

1.1 

5.93-10^ 

Chicago 

3.46-10-'^ 

1.1 

5.90-104 

San Francisco 

3.20-10-'^ 

1.02 

2.13-104 

Orlando 

3.64-10-'^ 

0.85 

6.31-10^ 


involve 114 cities and account for about 75 percent of all 48-state passengers and 70 percent of total domes¬ 
tic passengers. Data on the metropolitan population of the American cities are taken from the United States 
Census Bureau (2011 estimates). The travel rates are computed according to (4.1). 

Disease parameters within cities are taken from [18]. In that paper, the contact rate between susceptibles 
and infectives, j3*, is generally estimated to be 1.0. Like many respiratory diseases, influenza exhibits a 
seasonal pattern with a low summer and a high winter incidence. Here, we chose to only focus on the 
autumn-winter period going from October to March. Cities are divided into five general zones based on 
the average number of heating degree-days and cooling degree-days, see for instance [3]. We applied a 
scaling factor to j3* which depends on the climate zone of each city and ranges between 0.85 and 1.1 in 
the autumn-winter period. This provides a specific value j3, for each city i. The complete list of cities with 
their corresponding metropolitan population and transmission rate is provided in Table A.6 in Appendix F. 

Finally, the average length of the infectious period is estimated to be 2.95 days independently of the 
city [18], which leads to d, = d = I /2.95 for all i. 

4.3 Results and Discussion 

In this section, we use the results of Section 3 to compute the elasticity of the mean cumulative epidemic 
size D{t) and of the basic reproduction number Rq, with respect to the parameters of the model, over the 
first two weeks of the epidemic, that is, for f G [0,14] (as discussed in Appendix C, the branching process 
approximation is reasonable during this time period). 

We present our results for a small number of cities, namely New York, Orlando, Chicago, and San 
Francisco. As shown in Table 2, this set of cities covers all possible values of transmission rates. Table 2 
also presents the travel rates from New York to the three other cities, and the mean cumulative size of the 
epidemic after 14 days if it started in each of the four cities. We obtain f2o = 2.1401 and Rq = 3.2433, 
meaning that the epidemic breaks out with positive probability. Note that these values are particularly high 
due to the relatively high transmission rates during the autumn-winter period (between 0.85 and 1.1). 

Elasticity of the epidemic size with respect to the transmission rates. Figure 2 shows 
the elasticity of the mean cumulative size of the epidemic with respect to the transmission rate in the city 
where the disease started, for three cities with different transmission rates. We see that the elasticities are 
non-negligible and the curves appear in the same order as the transmission rates. Using formula (3.3), we 
see that an increase by 1% of the transmission rate in New York will induce an increase by approximately 
14% of the mean cumulative epidemic size after two weeks, while an increase by 1% of the transmission 
rate in Orlando will induce a corresponding approximate increase of 5% only. We observe that the elas¬ 
ticities increase over time except, interestingly, for Orlando where the elasticity reaches a maximum and 
then starts decreasing. This can be interpreted as an effect of the relatively low transmission rate associated 
with Orlando and the high travel rate from that city: if the epidemic starts in Orlando, after about ten days 
the growth of the disease in other cities connected to Orlando starts exceeding the infected population in 
Orlando (due to higher transmission rates in these cities, see also Figure C.l in Appendix C). There is thus 
a time from which the growth of the disease depends less on the transmission rate in Orlando than on the 
transmission rate in other cities, and therefore a small increase in the transmission rate of Orlando would 
start having a decreasing impact on the size of the epidemic. 
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Fig. 2. Elasticity of the mean cumulative epidemic 
size with respect to the transmission rates. Elastic¬ 
ities of the mean cumulative epidemic size at time t if 
the disease starts in city i (D/(r)) with respect to the 
transmission rate in city / (jS,) (left) for a sample of 
origin cities /. 



Fig. 3. Elasticity of the mean cumulative epidemic 
size with respect to the travel rates. Elasticities of 
the mean cumulative epidemic size at time t if the 
disease starts in city i (D/(0) with respect to the travel 
rate from city i to city j (Ti;), with New York as the 
origin city i of the disease, and for a sample of cities 

J- 


Table 3. Elasticity of the basic reproduction number Rq with respect to some model parameters which are constant for all cities. 

p: P* d c 

1.2357 -1.2350 -6.62-lO^^^ 


Elasticity of the epidemic size with respect to the travel rates. Figure 3 shows the elas¬ 
ticity of the mean cumulative number of infected individuals with respect to three travel rates from New 
York (where the epidemic is initiated): the travel rates to Chicago, to San Francisco, and to Orlando (see 
Table 2 for their values). We see that the elasticities are significantly smaller in absolute value than when 
they were computed with respect to the transmission rates. A small perturbation in the travel rates has 
therefore less impact on the dynamics of the disease than an equivalent perturbation in the transmission 
rates. We also observe that the elasticities are all negative, which suggests that increasing the travel rates 
out of New York, the origin city of the disease, may be benehcial from a sanitary point of view. This result 
might be surprising but it makes sense because the transmission rate in New York is the highest (j3, = 1.1), 
and we assumed that there is no transmission during travel (we refer to Appendix D.2 to see what happens 
when on-board transmission is taken into account). Indeed, increasing the travel rates from New York to a 
city with the same transmission rate, such as Chicago, will have almost no effect on the dynamics of the 
disease, while increasing the travel rates from New York to a city with a very low transmission rate, such 
as Orlando, will have a higher impact on the size of the epidemic. This argument further suggests that a 
smaller transmission rate in the destination city induces a larger elasticity (in absolute value), as confirmed 
in the figure. 

Elasticity oe Rq. Recall that the basic reproduction number Rq is an important threshold scalar quan¬ 
tity which measures the initial disease transmission on the whole network. In order to study the sensitivity 
of Rq with respect to the parameters of the model, we use the results of Section 3.4, where the expression 
for dpR (provided in (A.2) in Appendix A) simplifies to 

dpR = [l+T-^dmg{p)Y^ [T-\dpT)T-^dmg{p)-T-^dpdmg{^)] {[l + T-^dmg{p)y^ + 1 }. (4.2) 

Not surprisingly, the local sensitivity of Rq with respect to the transmission rate in one single city, j3,, 
or with respect to the travel rate between two specific cities, 7]j, is almost negligible (of the order of 10^® 
or less). We also performed a global sensitivity analysis of Rq with respect to some parameters which are 
assumed to be the same for all cities, namely the contact rate j3* = 1 and the removal rate d = 1/2.95. 
The results are shown in Table 3. We see that a 1% increase in the contact rate j3* (or, equivalently, in the 
transmission rate of all cities together) results in a 1.24% increase in Rq. We also note that, in absolute 
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value, the elasticities of Rq with respect to j3* and with respect to d are very close to each other, Rq being 
slightly more sensitive with respect to j3* than with respect to d. 

In order to analyze the relative change in Rq when all travel rates increase by the same factor, we write 
Tij = c ■ Tij (i 7 ^ j) where the constant c is equal to 1, and we compute the sensitivity of Rq with respect to 
c. The result is given in the last column of Table 3 and shows that if all travel rates increase by 1%, then 
Rq decreases by 6.62 • 10^^%. This suggests that increasing the frequency of all domestic travel (by the 
same percentage) would globally slow down the epidemic very slightly. This is again a surprising result 
which comes from the fact that each city has its own transmission rate, and increasing the frequency of all 
domestic travel (without on-board transmission) would, among other things, increase the travel rates out of 
the cities with a high transmission rate. Finally, we see again that the sensitivity is significantly lower with 
respect to the travel rates than with respect to the transmission rates. 

Remark 4.1 Note that the sensitivity results are very dependent on the disease under consideration: since 
the climate plays a fundamental role in the spread of influenza, the growth of the disease is sensitive to the 
travel rates because different cities in different climate zones exhibit different transmission rates. However, 
other types of infectious diseases may spread independently of the climate, in which case the air travel 
network would have much less impact on the epidemic growth, if we assume that there is no in-flight 
transmission. 

A discussion of the results in the case where on-board transmission is taken into account is given in 
Appendix D.2. 

4.4 Vaccination 

Vaccination can be used to control the spread of the disease within a population and hence eradicate it. 
From a public authority’s perspective, one may want to solve two questions: 

(i) what is the smallest proportion of individuals that need to be vaccinated in order to prevent an epi¬ 
demic outbreak, and 

(ii) once this number is reached, how many new infections do we prevent with each additional vaccinated 
individual? 

In this section, we will tackle those two issues by considering the simple problem where there is no on¬ 
board transmission and vaccination is uniform over the population. Appendices E.2 and E.3 respectively 
deal with the cases where vaccination is treated differently in each city and where there is on-board trans¬ 
mission. 

As we do not consider age patterns, we shall assume that each individual is vaccinated with some 
probability r, independently of his age and independently of the city. The parameter r thus also denotes the 
fraction of the total American population which is vaccinated. Vaccination reduces susceptibility, which, in 
our branching process approximation model, is equivalent to considering the new transmission rate vector 
(in the sequel, we write the subscript v each time a quantity depends on this new vector). 
The minimal fraction of individuals that should be vaccinated is solution of the minimization problem 

minimize rVv, 

r “ 

i 

subject to < 0, 

where Xmaxi^v) denotes the Perron-Erobenius (dominant) eigenvalue of the matrix The solution of 
this problem can be written out explicitly and the critical value obtained for our model is = 0.69, which 
means that at least 69% of the population has to be vaccinated to avoid an outbreak (see Appendix E.l for 
more details). Assuming that the total population in the US is 3.12 • 10*, the critical number of individuals 
to vaccinate is 3.12 • 10*rc « 2.16 • 10*; this solves Question (i). Given that at least that fraction of the 
population is vaccinated, the answer to Question (ii) is then given by the sensitivity of the mean total 
cumulative number of infected individuals until eradication with respect to r. 
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Table 4. Uniform vaccination case in the almost critical regime where r is such that 3.12-10^ (r — r^) = 10: the mean total cumulative 
size of the infection given the origin city of the disease (D/), and the number of infections prevented by introducing one additional 
vaccine (dDi), for three origin cities of the disease 

Origin city i Di dDi 

New York 1.1-lO^ 8.8-10^ 

San Francisco 2.6 • 10^ 2.0 • 10^ 

Orlando 3.3 lO^ 2.6 10^ 


Let Di represent the mean total cumulative size of the infection given that it was initiated by a first 
individual in city i, and D = (Di) = \\mt^ooD{t), where D{t) is defined in (2.10). Since the dominant 
eigenvalue of is negative, lim,^„oexp(f2vf) = 0, and from (2.10) we obtain 

D= (4.3) 

where D depends on r through f2y. Therefore, the sensitivity of D with respect to the fraction r of vacci¬ 
nated people is given by 


or 

= -{-^vy'diag{P){-Q^y'd. 

From this formula, we can approximate the mean number of prevented infections per additional vaccine 
given that a proportion r > rc of the population is vaccinated; one more vaccine leads to dr = (total 
population)^' = (1 /312) 10^®, so that 

do = -(1/3.12) 10^* (-f2y)-' diag(/3) d. 

Assume that we are in an almost critical regime, that is, the initial fraction r of vaccinated individuals 
is slightly larger than consider, for instance, a value r such that 3.12 • 10^ (r — r^) = 10, that is, such that 
the critical vaccinated population is increased by 10. The benefit of one additional vaccine in this regime is 
shown in Table 4, in which we compare the mean total cumulative size of the infection in the almost critical 
regime given the origin city i of the disease, Di, with the approximate number of prevented infections if 
we introduce one additional vaccine in the population, dDi, for three origin cities of the disease. We see 
that the introduction of one additional vaccine in this almost critical regime would reduce the size of the 
epidemic by a bit less than 10%. The number of prevented infections per additional vaccine is thus very 
large when the initial number of vaccinated people is close to the critical number; in other words, in the 
almost critical regime, the growth of the disease is highly sensitive with respect to the vaccination ratio. 
The values of dDi decrease very rapidly when the initial fraction of vaccinated people r becomes larger 
than the critical value r^. This is shown in Figure 4 for the three cities considered in Table 4, where we 
depict the number of people who would escape from the disease if we introduce one additional vaccine in 
the population as a function of the difference 3.12 • 10* (r — r^) between the initial number of vaccinated 
people and the critical number. 

5. Conclusion 

Explicit formulas are derived for the sensitivity analysis of a model of Markovian branching process evolv¬ 
ing on a fixed network. This tractable stochastic process is then used to model the early stages of a seasonal 
influenza-like disease speading on a network of cities in the United States. This approach allows us to study 
the sensitivity of the size of the epidemic and of the basic reproduction number with respect to the trans¬ 
mission and domestic travel rates in an analytic way, that is, without the use of simulations. Our analysis 
highlights the differences in the sensitivities with respect to the different parameters, confirming that a 
precise estimation of the transmission rates is far more important than a precise estimation of domestic 
mobility data. We also treat an extension of the epidemic model by considering vaccination campaigns. In 



SENSITIVITY ANALYSIS OF A BRANCHING PROCESS EVOLVING ON A NETWORK 


17 of 30 



Fig. 4. Uniform vaccination. Number of prevented infections per additional vaccine as a function of the difference between the 
initial number of vaccinated people (3.12 ■ 10** r) and the critical vaccinated population (3.12 ■ 10* r^), for three origin cities of the 
disease. 


this case, a sensitivity analysis enables us to calculate the marginal gain of one additional vaccine, show¬ 
ing how, at the almost critical regime, each additional vaccination prevents a large number of possible 
infections. 

Note that in the epidemic application, we chose to focus mainly on the sensitivity of expected quanti¬ 
ties, such as the mean epidemic size, which can also be described by a linearised SIR model, for simplicity 
and ease of comparison with those popular models. However our methodology also applies to more sophis¬ 
ticated quantities not deducible from a deterministic model, such as the variance of the epidemic size, or the 
probability that the epidemic dies out before a given time. We have considered the example of influenza, 
which is suitable for our purpose because people usually continue to travel when they are infectious, how¬ 
ever our model and methods are not restricted to this particular disease. 

Our methodology can be applied to a variety of other Markovian models involving reproduction, death 
and (discrete) node transition events. Note that the nodes of the network may not only represent geographic 
locations but also more abstract features such as physiological states of individuals, see for example [23]. 
A possible extension of the model would be to generalize the results of Sections 2 and 3 to a dynamic 
network (instead of a fixed one); this is the topic of ongoing research. 
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A. Further results on the branching process on a network and its sensitivity 

In this section we gather some technical results related to Sections 2 and 3. 

A. 1 The matrix R 

Recall from Section 2.4 that the entry of the matrix R is the expected total number of children born at 
node j from a parent born at node i, during the entire lifetime of the parent. For each k ^ 1, define the 
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n X matrix 

% = -T-'Bk, 

whose entry {^k)i;jij 2 ---jkJ records the probability that an individual at node i eventually gives birth before 
dying (potentially after moving to other nodes), and the first reproduction event involves the birth of k chil¬ 
dren and the parent moving to node j, while the children start their life at nodes ji, 72 , ■ ■ •, jk respectively. 
Note that nodes j, 71 , 72 , ■ ■ • j jk are not necessarily adjacent to node i here. An explicit expression for R can 
be given in terms of the matrices as shown in the next proposition. 


Proposition A.1 The matrix R of expected total progeny size is given by 


R = 




k-l 


i=0 


(A.l) 


Proof. We develop the entry Rij by conditioning on the first reproduction event happening to the parent at 
node i: 


k^l h Jk ( 

where 5 ,7 = 1 if and only if i = j. In matrix form, this becomes 




k>l 


k-l 


E (1^'^ 0/(g) +(1^0/)/? 


1=0 


and since the matrix Y.k^i (g)/) is substochastic, we can write R explicitly as in (A.l). □ 

Using the explicit expression for R given in (A.l), the sensitivity of R with respect to a parameter p is 
given by 

^ k^l ^ ^ k^l k^l i=0 

(A.2) 

k^l k^l 1=0 

where dp% = dp[-T-^Bk] = -T-^[dpT % + dpBk]. 


A.2 The sensitivity of Q{t,s) and G{t,s) 

Here we are interested in characterizing dpQ{t,s) and dpG{t,s). By differentiating (2.5) with respect to p, 
we obtain 


where 


d,dpQ{t,s) = dpd + dpTQ{t,s) + TdpQ{t,s)+Y,dpBkQ{t,s)^’^+^ 

k 

k i=0 

= A(t,s)+B(t,s)dpQ(t,s) 

A(t,s) = dpd + dpTQ(t,s)+J^dpBki2(f,s)^^+^^ 


B(t,s) = T + J^Bk£(Q(t,s)(‘)0/0Q(t,sf-‘^). 

k 1=0 


(A.3) 


Unfortunately the matrix linear differential equation (A.3) for dpQ{t,s) does not have any closed-form 
solution because the matrices B{t,s) do not necessarily commute for different values of t. 

Using the same argument and Eq. (2.9), we show that dpG{t,s) satisfies the same matrix linear differ¬ 
ential equation as dpQ{t,s), the only difference being that the term dpd is replaced by dpds. 
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A.3 An approximation for the sensitivity ofM{t) 

An alternative method to compute dpM{t) provides a new approximation for dpM{t) which relies on the 
integral representation of the derivative of the matrix exponential [34]: 

dpM{t) = dpe\p{Qt) = [ exp(f2T)(9„f2 exp(f2(f — t)) t/T. (A.4) 

Jo 

Proposition A.2 The matrix dpM{t) can be approximated as 

1 K(t)K(t) + \ (ct)^ 

- exp(f2of)Av ^ Y, ew{-ct)^-rl-&‘Af^dpnAv©’'^‘^^Af^ (A.5) 

^ i=0 k=i+l 

where c, Ay, K{t) and 0 are defined in the proof. 

Proof. The integral in (A.4) is approximated by using a duality argument (as in [24]) and the uniformiza- 
tion technique for continuous-time Markov chains (see for instance [36]), as we detail now. 

Let and v be the left and right eigenvectors of Q. corresponding to the Perron-Frobenius eigenvalue 
f2o, normalized by v^l = 1 and u^v = 1, and let Ay = diag(v). We first rewrite exp(f2f) as 

exp(f2r) = exp(f2or)Av exp(0r)Ay^*, 


where 

It is easy to show that 0 satisfies all the properties to be the generator of a continuous-time Markov chain; 
©ij 0 for i f i and 0 ©a = — - This Markov chain is called the dual of the branching process 

with mean population size matrix exp(f2t). The integral in (A.4) then becomes 

J exp{QT)dp£2 exp{Q{t — T))dT = 

exp(f2ot)Av / exp{©z)Af^dp£2Avexp{©{t — T))dxAf^. 

Jo 

Since exp(0f) is now a probability transition matrix, we can use the uniformization method to solve the 
integral; let c > max, |0„| and 0 — I +^. We can write 

exp(0f)=I:exp(-cf)^0^ 

and since the matrix 0 is substochastic, we have 

||exp(0O|| ^ i:exp(-cO^||0'=|| 

k^O 

/ , , (ct)^ 

k^O 

Let K{t) be such that Lfio exp(—cf) ^ 1 — e for a given e > 0. Then, 

K{t) (pf\k 

exp(0f) « E exp(-cf) 0^ 
k=0 


and by replacing exp(0T) and exp(0(f — t)) in the last integral by their approximations, we finally obtain 
the approximation (A.5) for dpM{t). □ 
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Fig. B.l. Elasticity. Interpretation of the elasticity in terms of percent of increase or decrease of a measure of the model 


B. Interpretation of the elasticity 

In this section we interpret the elasticity in terms of a proportional response to a proportional pertur¬ 
bation. Let a be the elasticity of a measure X of the model with respect to a parameter p, that is, 
Tp = d\ogX/d\Qgp = a. Assume that p increases by r% (r relatively small), and let z be the induced 
proportional increase (or decrease, depending on the sign of a) of X, that is, 

p' = p(\ -fO.Olr) =X(1 +z). 


By taking the logarithm we obtain 

logp' = logp + log(l -l-O.Olr) \ogX' = logX + log(l +z). 


We thus have 

_ 5logX ^ log(l+z) 
d\ogp log(l+0.01r) ’ 

and we obtain 

z « exp(fllog(l +0.01r)) — 1. 

So, for instance, if p increases by r = 1 %, then thenX increases (or decreases) by approximately 100[exp(fl/100) — 
1]%. In Figure B.l, we show the curve y = 100[exp(a/100) — 1] and we compare it with y = a. We see 
that the two curves are almost superimposed when |fl| is small. 


C. Validity of the branching process approximation 

In this section, we compare the mean instantaneous and cumulative population size obtained from the 
branching process defined in Section 2 to their analogue in the popular deterministic SIR model. 

For that purpose, we illustrate the two measures M{t) and D{t) on an example where there is a first 
infectious case in Orlando (city i) at time f = 0. This city has a low transmission rate (j3, = 0.85) but it 
is well connected with New York (city j, Tij = 0.0033) which has a high transmission rate (jSy = 1.1). 
In Figure C.l, we show the growth of the population of infected individuals within the origin city of the 
disease, as well as in New York. In that figure, we compare the branching process approximation to the SIR 
model {{Si{t),Ii{t),Ri{t)), 1 ^ ^ 114} which tracks, for every city i, the number of susceptible (5,(f)), 
infected and removed (R,(f)) individuals in city i at time f, and whose evolution is modeled by the 
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Fig. C.l. Mean instantaneous and cumulative epidemic sizes. Mean epidemic sizes between the first and the 35th day (in log¬ 
arithmic scale) if the disease is initiated by one infectious case in Orlando (city i) on day t — 0. The plain lines correspond to the 
branching process approximation, the dashed lines correspond to the SIR model. M„ (0 is the mean instantaneous number of infected 
individuals in Orlando, Mij{t) is the mean instantaneous number of infected individuals in New York (city j). Mi, {t) — is the 

mean instantaneous total number of infected individuals in all cities, and A (0 is the mean cumulative number of infected individuals 
in all cities. The branching process approximates the SIR model reasonably during the first three weeks of the epidemic. 


Mi Mi 

TT Si U + L ^7 Tji — li 'Y^Tij — di li 
Mi Mi 

diU + L^ i'^ji ~ ^i V Tjj, 

Mi Mi 

where {Si{Q),Ii{Q),Ri{0)) is specified, for 1 ^ ^ 114. 

We assume homogeneous mixing within each city. According to [11], during the course of a major 
epidemic, the disease grows like a branching process within each city i until about ^/Ni members of the 
population of A, individuals become infected. If we use this criterion to determine the maximum validity 
period of the branching process approximation in each city, we find that the approximation is accurate 
during the first 10 days of the epidemic on average, the time value ranging from 6 days (Aspen, Colorado, 
VA; = 81.6, Pi = 1.1) to 15 days (Houston, Texas, VAj = 2467.2, pi = 0.85). 

This criterion is probably too strong for our purpose. If we rather fix to 1% the relative error obtained 
when computing with the branching process versus Ii{t) with the SIR model, the approximation is 

accurate until 13 days on average. The maximum approximation time depends on the metropolitan size 
of the city considered and on its associated transmission rate; it ranges from 5 days (Aspen) to 20 days 
(Dallas, Texas, ^/Ni = 2554.8, pi — 0.85). For Orlando (VA, = 1473.4) and New York (\/Ni = 4366.9), 
the two cities used to illustrate our measures, the approximation is accurate until respectively 19 and 15 
days after the introduction of the first infectious case in the city. 

Figure C.l also shows that, while the epidemic initially started in Orlando, after a bit more than two 
weeks it mainly evolves in New York because of the relatively high travel rate between Orlando and New 
York and the high transmission rate in New York. This is in line with [17] and confirms that air travel has 
a non-negligible effect on the spatial spread of the disease. 

Remark C.1 The expected instantaneous epidemic sizes Mij{t) obtained from the branching process 
approximation correspond to the solution of the linearisation of the quadratic differential equation for Ij{t) 
around the disease-free equilibrium, that is, if we assume that Sj « Ay for all j. Note that, thanks to its 
stochastic nature, the branching process approximation may also allow us to study other types of questions. 


nonlinear differential equations 

dSi 

dt 

dli 

dt 

dRj 

dt 
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such as the full distribution of the instantaneous and cumulative epidemic sizes during the early stages of 
the disease. These questions would be particularly relevant if the epidemic starts with a small number of 
infected individuals, as the strength of stochasticity increases as population sizes get smaller. 


D. On-board transmission 

In this section, we show what the results of Sections 4.1 and 4.3 become when on-board transmission is 
taken into account. 

D.l Description of the model 

To the best of our knowledge, there is not enough available data to conclude in an acceptable on-board 
transmission rate for influenza-like diseases. Therefore, for the sake of simplicity, we assume that during 
a flight an infected individual may infect up to two susceptible individuals (his/her two closest neighbours 
in the airplane), independently of each other. This is probably a lower bound of the real number of new 
infections, but it is sufficient for our purpose of showing the differences with the case without on-board 
transmission. 

We model the number of new infections generated by an infected individual during a flight from city i to 
city j with a binomial distribution B(2,py ), where ptj is the probability that an infected individual infects 
another individual during a flight from city i to city j. Since no real data are available for these on-board 
transmission probabilities, for the sake of our analysis we shall assume that they take the hypothetical form 

5hij 

5hij + V 

where hij is the flight time (in hours) from city i to city j, which is directly proportional to the distance 
between i and j. The data on non-stop distances between cities are available from the US Department of 
Transportation [2]. For i f j, the probability {Pi)ij that an infected individual infects exactly one suscep¬ 
tible individual during the flight between city i and city j is given by {Pi)ij = 2p,;(l — Pij)- Similarly, the 
probability {P 2 )ij that an infected individual infects exactly two susceptible individuals during the flight 
between city i and city j is given by {P 2 )ij = pjj- These probabilities depend on the flight time and are 
shown in Figure D.l. 

The travel rate from city i to city j without subsequent transmission, denoted by Tip or associated with 
the transmission of the disease to one or two individuals, denoted by {Cfjij and (C 2 ),; respectively, are 
given by 

Th = Tp[\-{P,+P2)ij\, 

{Cx)ij = 7], (Pi)/;, 

{C2)ij = Tp{P2)ij, 

for i f j. 

In the branching process approximation, the travel rate matrix then becomes T' (the diagonal of T' 
being the same as that of P), and if we assume that air travel is instantaneous (that is, flight times are 
negligible), there may be up to two transmissions at a time. Therefore, the birth rate matrices are such that 
the only nonzero entries of the matrix B\ are 

(P'l )/;,; = A- and (p'i),..,.,. = (Ci),;, 


the only nonzero entries of the matrix Pj ^6 


= (C2) 




and 


BI = 0 for k ^ 3. 



SENSITIVITY ANALYSIS OF A BRANCHING PROCESS EVOLVING ON A NETWORK 


23 of 30 



Fig. D.l. On-board transmission probabilities. Probability that an infected individual infects one or two new individual(s) as a 
function of the flight time (in hours). 


Table A.5. Mean cumulative epidemic sizes after 14 days associated with four origin cities if on-board transmission is taken into 
account. We indicate the ratio of the mean epidemic size with on-board transmission to the corresponding mean size without on¬ 
board transmission as given in Table 2. 


City j 

D,(14) 

Ratio 

New York 

6.49 • 10^ 

1.0945 

Chicago 

6.70-10'^ 

1.1347 

San Francisco 

3.25-10'^ 

1.5231 

Orlando 

1.71-10'^ 

2.7003 


With these, Bj (Z® 1) = Bj(l 0/) = diag(j3) +Ci and ^ 2 . 

so that the matrices Q. and R become 

Q' = 7’'+2(diag(/3)+Ci) + 3C2 

R’ = -[/+r'-'(diag(/3)+Ci+C2)]^'r'-'(diag(/3)+Ci+2C2). 


Note that our analysis could be adapted to more than two new infections during a flight. This would 
require to define subsequent matrices Pi and C, for i > 2, but would complicate the model unnecessarily for 
our purpose. 

D.2 Results and Discussion 

We found that the difference in the global growth of the disease (that is, in Rq) is very small when taking 
on-board transmission into account. However, at the city level, the difference can be more noticeable, as 
shown in Table A.5: for example, we see that if the disease starts in Orlando, the expected cumulative 
epidemic size after two weeks is more than twice larger with on-board transmission than without. 

In order to emphasize the effects of the on-board transmission mechanism on the elasticities, we shall 
compare elasticities without and with on-board transmission on the same graphs, even if for some cities the 
branching process approximation might be less accurate around f = 14 days when on-board transmission 
is taken into account. The plain lines in the graphs correspond to the model with on-board transmission, 
and the dashed lines correspond to the case without on-board transmission. 


Elasticity of the epidemic size with respect to the transmission rates. In Figure D.2 
we show a comparison of the elasticities of the mean cumulative size of the epidemic with respect to the 
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Fig. D.2. Elasticity of the mean cumulative epi¬ 
demic size with respect to the transmission rates. 

Same graphs as in Figure 2, with on-board transmis¬ 
sion; the plain lines correspond to the model with on¬ 
board transmission, and the dashed lines correspond 
to the model without on-board transmission. 



Fig. D.3. Elasticity of the mean cumulative epi¬ 
demic size with respect to the travel rates. Same 
graphs as in Figure 3, with on-board transmission; 
the plain lines correspond to the model with on-board 
transmission, and the dashed lines correspond to the 
model without on-board transmission. 


transmission rates with and without on-board transmission. We see that in general on-board transmission 
decreases markedly the elasticity of £>, (?) with respect to j3, . 

Elasticity of the epidemic size with respect to the travel rates. In Figure D.3, we see 
that the elasticity of the mean cumulative size of the epidemic with respect to the travel rates is larger 
when on-board transmission is taken into account. This indicates, as expected, that in that case a change 
in the travel rates would affect much more the dynamics of the disease than when on-board transmission is 
not taken into account. More precisely, we see that the elasticities start increasing to positive values with 
respect to the travel rates out of New York, due to the risk of transmission on board airplanes. After a 
few days, however, the elasticities reach a maximum value and start decreasing to negative values, which 
is especially clear for Orlando. This indicates that, in the long term and under our on-board transmission 
assumptions, it would still be beneficial from a sanitary point of view that infected people travel out of New 
York to cities where the transmission rate is lower. 

Elasticity oe Rq. Let C = Ci -I- C 2 . The expression given in (A.2) for dpR simplifies to 

dpR' = [i+r-\dmg{p)+c)Y' [r-\dpr)r-\dmg{p)+c)-r-'dp{dmg{p)+c)] ■ 

{[/ + r'-'(diag(j3)+C)]^'+/}. (D.l) 

The values of the elasticity of Rq with respect to the different parameters with on-board transmission are 
identical to those obtained in the case without on-board transmission up to the 6th decimal. We conclude 
that on-board transmision (as we define it) has a negligible effect on the sensitivity of Rq with respect to 
these parameters. 


E. Vaccination 

The vaccination problem is stated as follows: given a population and an epidemiological model, what is 
the smallest fraction of the population that needs to be vaccinated to prevent the epidemic to break out? 
And once this number is reached, how many new infections do we prevent with each additional vaccinated 
individual? We first answer these two questions in the scalar case (uniform vaccination) and the vector case 
(city-dependent vaccination) without considering on-board transmission. We then discuss what the results 
become when on-board transmission is taken into account. 


E.l Uniform vaccination 

Minimal fraction to be vaccinated. In our case, if we assume that vaccination is done uniformly 
over the population of each city, without consideration of age, social status, etc., then vaccination has the 
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same effect as reducing the transmission rate. Let us suppose for now that the vaccination is done in the 
same proportion over the whole country, and that a fraction r of the population is vaccinated. It leads to a 
new transmission rate vector P, = il-r)P. In the sequel, we write the subscript v each time a quantity 
depends on this new vector. The minimal fraction of the population to be vaccinated in order to prevent a 
breakout of the epidemic is then given by the minimization problem: 


minimize r V Ni 

r “ 

i 

subject to Xmaxi^v) ^ 0. 


(E.l) 


If we write, for simplicity, = diag(/3), one has Q = T + 2Ap and Qy = T + {2 — r)Ap. Since 
Qy is not symmetric (it results from the sum of T, which is not symmetric, and a diagonal correction 
{2 — r)Ap), there is no guarantee that its eigenvalues are real. However, since T = N^^A, where N is 
diag(A^,) the diagonal matrix of metropolitan populations, and A is a symmetric air travel matrix (as in 
[17]), the matrix is also symmetric and has the same eigenvalues as T. Consequently, Qy 

has the same eigenvalues as N^I^QyN^^/'^ which is symmetric too. For symmetric matrices, the condition 
Xmax{^\) ^ 0 can be rewritten using Rayleigh quotients as 


v^N^I^QyN-'^l^v 

yTy 


^ 0 , 


Vv : ||v||^ > 0. 


By expanding the expression of Qy and simplifying, one obtains 

{r-2)v^Apv, Vv : Hvf > 0. 

— 1/2 

If all j3, are nonzero, then Ap is invertible, and we can write v — Ap ' vv, so that the condition becomes 


w'^Ap ^/^N'^I^TN-'-I^Ap < (r-2)w"^w, Vw : ||wf > 0. 

By dividing both sides by w^w, one obtains 

- ^ -^(t-2), Vw:|1m'|/>0, 


which can finally be translated into 

2 +Xmax{Ap^T) ^ r. (E.2) 

In our case, we conclude that r should be larger than 0.6917, so at least 69.17% of the population has to be 
vaccinated in order to prevent the epidemic to break out. 


E.2 City-dependent vaccination 

Minimal fraction to be vaccinated. When vaccination is not done uniformly over the population 
but varies from city to city, the vaccination ratio is described by a vector r, where r, is the fraction of the 
population of city i that is vaccinated. If we define V = diag(r), the optimisation problem becomes 


minimize ^ r,A( 

^ i 

subject to Xmaxi^v) ^ 0, 


(E.3) 


with Qy = T (21 — V)Ap. Again, this matrix is not symmetric but it has the same eigenvalues as 
which is symmetric. The constraint Xmax(^v) ^ 0 can then be written as 


yT^l/2(7’ + (2/_y)4^)At-V2^,<^0 Vv: |lv|/ >0. 
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By writing again w = A 


1/2 

P 


V, one obtains after some elementary manipulations 


U,T[y _^1/2 m + 


>0 Vw ; IIm'Ip > 0, 


(E.4) 


which simply means that V — N^^^{Ap ^^'^TAp + 2/)A^ be positive semi-definite. 

In general, the solution of such an optimization problem is not trivial, but it can be solved using Linear 
Matrix Inequality (LMI) solvers. However, in this case, we can derive a fairly accurate approximation. 
Recall that the matrix T contains on its off-diagonal elements the travel rates between cities, while the 
diagonal elements are given by 7/,- = — jS,- — dt— ^ 7]^. In the present situation, the travel rates are small 

j^i 


with respect to transmission and recovery rates, and we can approximate T ss diag(—j3 — d). Then, the 
condition becomes that V — Ap^T — 27 is positive semi-definite. Since we want to minimize a weighted 

sum of the elements of V, with only positive weights, this reduces to solve r, — (—j3, — t/,) — 2 ^ 0 for 

each i, which immediately gives 


G > 1 - {di/^i). 


(E.5) 


We solved the original problem with condition (E.4), using the Matlab software CVX [19,20] for dealing 
with the LMI condition, and we denote the solution by rcvx- The values of the vector rcm range between 
0.5965 and 0.6929, with a mean of 0.6574. We compared this solution with the solution obtained in (E.5), 
denoted by rappmx- In all cities the optimal vaccination fractions obtained with both methods are very 
similar: the largest relative error between the entries of Vcvx and rappmx is 1%. There is a clear correlation 
between the intensity of the travel rates and the relative error; the approximation is excellent in weakly 
connected cities and less accurate in the most connected cities. 

However, this approximation is valid in our case because travel rates are very small with respect to the 
transmission and removal rates J3 and d. This is no longer the case when travel has a larger impact than 
here, that is, when commuting traffic is added, or if on-board transmission is taken into account, as shown 
in Appendix E.3. 


Sensitivity of D with respect to r. Now, the sensitivity of vector D with respect to r, (the other 
rj for j 7 ^ i being held constant) is given by 


dD 

dvi 




(9f2v 

dr-i 




where the only nonzero entry of the matrix dLl^ldri is the entry (i, i) equal to —j3,, and dri = 1 /A,. In this 
case, when the proportion r of vaccinated individuals in each city is large enough to prevent an epidemic 
to outbreak, the effect of an additional vaccine in city i is then given by 


dD=—(-Qx 

Ni 




Similar to the uniform case, the value of dDi decreases rapidly when r moves away from r^, as shown 
in Eigure E.l, which depicts the number of people who would escape from the disease if we introduce 
one additional vaccine in the origin city of the disease, as a function of (r — Tc); x Ni, the initial additional 
number of vaccinated people with respect to the critical number in the city, (rc),A,. 

Remark E. 1 If we assume that the vaccination campaign only applies to the 114 cities considered, then the 
results in this section show that, in order to prevent the epidemic to break out, at least g I], A/ = 1.64 • 10* 
individuals have to be vaccinated in the uniform case, as opposed to (fe), A, = 1.57 • 10* individuals in the 
city-dependent case. This highlights the benefit of doing city-dependent vaccination, since approximately 
seven million vaccinations less are needed to avoid an epidemic outbreak. 


E.3 Effect of on-board transmission 

The optimal vaccination fraction can be estimated in a similar way when on-board transmission is taken 
into account. The optimization problems (E.l) and (E.3) stay unchanged, except for the expression of Ely, 
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Fig. E.l. City-dependent vaccination. Number of prevented infections per additional vaccine in the origin city of the disease as a 
function of the difference between the initial number of vaccinated people in the city (r, V/) and the critical vaccinated population in 
the city ((rc)/V/), for three origin cities of the disease. 


which is now given by 

i2, = r' + (2/-y)4^+2Ci+3C2. 

For simplification reasons, we write T\ —T'-\- 2C\ + 3C2, so that = 7i + (2/ — V)Ap. 

Note that similar to T, T\ is not symmetric, but can be written as the product of and a symmetric 
matrix, so that T\ has the same eigenvalues as which is a symmetric matrix (the same holding 

for r'. Cl and C 2 ). 

Uniform vaccination. The same manipulation as in the case without on-board transmission provides 
the condition equivalent to (E.2): Xmaxi^v) ^ 0 if and only if 


2 + Xmax (4^ ' Ti) < r. (E.6) 

In our case, r should be larger than 0.6960, which is a bit larger than the minimal vaccination fraction 
without on-board transmission; this shows again the influence of on-board transmissions on the size of the 
epidemic. 

City-dependent vaccination. Similar to the city-dependent case without on-board transmission, the 
condition Xmaxi^v) ^ 0 of the optimisation problem (E.3) can be finally written as that 

y _ At 1 /2 ( 4 ^ 1/2 1 /2 2/)A?-1 /2 


has to be positive semi-definite. 

The optimisation problem can be solved using CVX and be compared with the simplified solution r approx 
obtained when assuming that the travel rates are negligible. The obtained vaccination rates Tcvx are higher 
than in the case without on-board transmission: their values range between 0.6015 and 0.7824 with a 
mean of 0.6682. Strikingly, since the influence of air travel is enhanced by the on-board transmissions, we 
observe larger differences between the solution rcvx returned by CVX and the simplified solution Vapprox- 
In the present case, 55 cities have a relative error between Tcvx and rapprox larger than 1% (the largest 
relative error is 12%), and there is a very clear trend between the connectivity of a city and the error; the 
most connected cities have a large relative error, while the most isolated cities have a well-approximated 
vaccination rate. 
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F. Complete data 

Table A.6 provides the complete list of the 114 American cities considered in this paper, with their corre¬ 
sponding metropolitan population (2011 estimates of the United States Census Bureau), and automn-winter 
transmission rate. The four cities which illustrate our sensitivity analysis in Section 4.3 are in bold. 

The 114 X 114 travel rate matrix is too large to be represented as a whole and can be computed according 
to (4.1) using the average daily number of passengers for each city pair obtained from the Domestic Airline 
Fares Consumer Report of the US Department of Transportation for the first Quarter of 2011 [2]. 
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Table A. 6 . The 114 cities considered, together with their metropolitan population Ni (based on the 2011 
estimates of the United States Census Bureau), and their winter transmission rate j3, (based on [17]). 


City i 

Ni 

A 

City i 

Ni 

A 

Albany 

857592 

1.1 

Little Rock 

685488 

1.02 

Albuquerque 

857903 

1.02 

Los Angeles Metro Area 4/ 

12870000 

1.02 

Allentown/Bethlehem 

816012 

1.1 

Louisville 

1259000 

1.02 

Amarillo 

246474 

0.85 

Lubbock 

276659 

0.85 

Aspen; CO (urban) 

6658 

1.1 

Madison 

570025 

1.1 

Atlanta 

5475000 

1.02 

Medford 

201286 

1.02 

Atlantic City 

271712 

1.02 

Memphis 

1305000 

1.02 

Austin 

1705000 

0.85 

Miami/Ft. Lauderdale 4/ 

5547000 

0.85 

Baltimore/Washington 4/ 

2691000 

1.02 

Midland/Odessa 

266941 

0.85 

Baton Rouge; LA 

786947 

0.85 

Milwaukee 

1560000 

1.1 

Bellingham 

200434 

1.1 

Minneapolis 

3270000 

1.1 

Billings; MT 

154553 

1.1 

Mission/McAllen/Edinburg 

741152 

0.85 

Birmingham 

1131000 

1.02 

Moline 

379066 

1.1 

Bloomington/Normal 

167699 

1.1 

Myrtle Beach 

263868 

1.02 

Boise 

606376 

1.1 

Nashville 

1582000 

1.02 

Boston/Providence 4/ 

6190000 

1.1 

New Orleans 

1190000 

0.85 

Buffalo 

1124000 

1.1 

New York Metro Area 4/ 

19070000 

1.1 

Burlington 

208055 

1.1 

Newburgh/Poughkeepsie 

677094 

1.1 

Cedar Rapids/Iowa City; lA 

450462 

1.1 

Newport News/Williamsburg 

1674000 

1.02 

Charleston 

659191 

0.85 

Norfolk 

1647000 

1.02 

Charlotte 

1746000 

1.02 

Oklahoma City 

1227000 

1.02 

Chicago Metro Area 4/ 

9581000 

1.1 

Omaha 

849517 

1.1 

Cincinnati; KY 

2172000 

1.02 

Orlando 

2082000 

0.85 

Cleveland/Akron 4/ 

2790935 

1.1 

Palm Springs; CA 

4143000 

0.85 

Colorado Springs 

626227 

1.1 

Panama City; FL 

164767 

0.85 

Columbia; SC 

744730 

0.85 

Pasco/Kennewick/Richland; WA 

245649 

1.1 

Columbus 

1802000 

1.1 

Pensacola; FL 

455102 

0.85 

Corpus Christ! 

416095 

0.85 

Philadelphia 

5968000 

1.1 

Dallas/Fort Worth 4/ 

6448000 

0.85 

Phoenix 

4364000 

0.85 

Dayton 

835063 

1.1 

Pittsburgh 

2355000 

1.1 

Denver 

2552000 

1.1 

Plattsburgh; NY 

81618 

1.1 

Des Moines 

562906 

1.1 

Portland 

2242000 

1.02 

Detroit 

4403000 

1.1 

Raleigh/Durham 

1627228 

1.02 

Eagle; CO 

61699 

1.1 

Reno 

419261 

1.1 

El Paso 

751296 

0.85 

Richmond 

1238000 

1.02 

Eugene; OR 

351109 

1.02 

Rochester 

1036000 

1.1 

Fargo; ND 

200102 

1.1 

Sacramento 

2127000 

1.02 

Fayetteville; AR 

464623 

0.85 

Salt Lake City 

1130000 

1.1 

Flint 

424043 

1.1 

San Antonio 

2072000 

0.85 

Fort Myers 

586908 

0.85 

San Diego 

3054000 

1.02 

Fresno; CA 

915267 

0.85 

San Francisco/Oakland 4/ 

4318000 

1.02 

Grand Rapids 

778009 

1.1 

Santa Barbara; CA 

407057 

0.85 

Greensboro/High Point 

714765 

1.02 

Santa Rosa; CA 

472102 

0.85 

Harlingen/San Benito 

396371 

0.85 

Sarasota/Bradenton 

688126 

0.85 

Harrisburg 

536919 

1.1 

Savannah; GA 

343092 

0.85 

Hartford 

1196000 

1.1 

Seattle 

3408000 

1.1 

Houston 

5867000 

0.85 

Sioux Falls; SD 

238122 

1.1 

Huntsville 

406316 

1.02 

Spokane 

468684 

1.1 

Indianapolis 

1744000 

1.1 

St. Louis 

2829000 

1.02 

Islip 

19070000 

1.02 

Syracuse 

646084 

1.1 

Jackson/Vicksburg 

589041 

1.02 

Tallahassee; FL 

360013 

0.85 

Jacksonville 

1328000 

0.85 

Tampa 

2747000 

0.85 

Kansas City 

2068000 

1.02 

Tucson 

1020000 

1.02 

Key West; FL 

73165 

0.85 

Tulsa 

929015 

1.02 

Knoxville 

699247 

1.02 

West Palm Beach/Palm Beach 

5547000 

0.85 

Las Vegas 

1903000 

0.85 

White Plains 

19070000 

1.1 

Lexington 

470849 

1.02 

Wichita 

612683 

1.02 



